THE CALCULATION OF THE DISTANCE TO A NEARBY 
DEFECTIVE MATRIX 



MELINA A. FREITAG* AND ALASTAIR SPENCEt 

Abstract. In this paper a new fast algorithm for the computation of the distance of a matrix to 
a nearby defective matrix is presented. The problem is formulated following Alam & Bora (Linear 
, Algebra Appl., 396 (2005), pp. 273-301) and reduces to finding when a parameter-dependent matrix is 

singular subject to a constraint. The solution is achieved by an extension of the Implicit Determinant 
Method introduced by Spence & Poulton (J. Comput. Phys., 204 (2005), pp. 65-81). Numerical 
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results for several examples illustrate the performance of the algorithm 
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1. Introduction. Let A be a complex n x n matrix with n distinct eigenvalues. 
It is a classic problem in numerical linear algebra to find 



d(A) — inf { || ^4 — i?||, B is a defective matrix}, 

where || • || = || • \\p or | • || = || • I2. Hence d(A) is the distance of the matrix A to 
the set of matrices which have a Jordan block of at least dimension 2. In this paper 
we present a fast numerical method to find a nearby defective matrix. We formulate 
the problem as a real 3-dimensional nonlinear system which is solved by Newton's 
method. Though not guaranteed to find the nearest defective matrix, since Newton's 
method provides no such guarantees, in all the examples considered our method did, 
1 in fact, find the nearest defective matrix and hence d(A) was computed. 

The distance of a matrix to a defective matrix is linked with the sensitivity analysis 
of eigenvalues. The condition number of a simple eigenvalue A is given by l/\y H x\, 
(see [12] ) where x and y are normalised right and left eigenvectors corresponding to 
A. For a defective eigenvalue we have y H x = and hence the condition number of 
the eigenvalue is infinite. But it is well-known that even if the eigenvalues of a matrix 
are simple and well-separated from each other, they can be ill-conditioned. Hence the 
measure of the distance d(A) of a matrix A to a defective matrix B is important for 
determining the sensitivity of the eigendecomposition. There is a very informative 
discussion and history of this problem in [2] , where the contributions of Demmel [3l [4] 
and Wilkinson [13j [14] are discussed in detail. Another important paper is that by 
Lippcrt and Edelman [8] who use ideas from differential geometry and singularity 
theory to discuss the sensitivity of double eigenvalues. In particular, they present a 
condition that measures the ill-conditioning of a matrix with a 2-dimensional Jordan 
block. The key paper that provides the solution to the nearest defective matrix 
is that of Alam and Bora [1] who provide both theory and an algorithm based on 
pseudospectra. 

Following Trefethen and Embree [llj , the e-pseudospectrum A e (A) of a matrix A 
is given by 

A E (A) = {a min (A~zI) < e}, 
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where e > and <J m i n denotes the smallest singular value. Equivalently 

A e (A) = {z e C | det(A + E - zl) = 0, for some E e C nxn with ||J5|| < e}. 

If A £ (A) has n components, then A + E has n distinct eigenvalues for all perturbation 
matrices E G C nx ™ with \\E\\ < e and hence A + E is not defective. Alam and 
Bora [T] take these ideas and seek the smallest perturbation matrix E such that the 
pseudospectra oi A + E coalesce. They present the following theorem (see [I] Theorem 
4.1] and [2J Lemma 1]). 

Theorem 1.1. Let A e C nxn and z e C \ A (A), so that A - zl has a simple 
smallest singular value e > with corresponding left and right singular vectors u and 
v such that (A — zl)v — eu. Then z is an eigenvalue of B = A — euv 11 with geomet- 
ric multiplicity 1 and corresponding left and right eigenvectors u and v respectively. 
Furthermore, if — 0, then z has algebraic multiplicity greater than one, hence it 
is a nonderogatory defective eigenvalue of B and \\A — B\\ = e. 

Proof. As A — zl has unique smallest singular value e > with corresponding left 
and right singular vectors u and v, we have rank(A — zl — euv H ) = n — 1 and hence 
z is an eigenvalue of B = A — euv H with geometric multiplicity 1. Further 

Bv = zv and u B = zu H , 

and u and v are left and right eigenvectors such that u H v = 0. So z is a multiple 
eigenvalue of B and hence z is a nonderogatory defective eigenvalue of B. □ 

Theorem 11.11 leads to the result E := —euv H so that B = A + E is a defective 
matrix and d(A) = e, since v H v — u u — 1. One drawback of the algorithm in [I] is 
that it is rather expensive since it involves repeated calculation of pseudospectra. Also 
a decision on when two pseudospectral curves coalesce is required. In [5] a method 
based on calculating lowest generalised saddle points of singular values is described. 
This has the advantage that it is able to deal with the nongeneric case when A — euv H 
is ill-conditioned. We shall present a straightforward, yet elegant and very fast method 
that deals with the generic case when A — euv H is well-conditioned. 

Using the notation of Theorem 11.11 the problem is to find z e C, u, v £ C" and 
e£l such that 

(A- zI)v-eu = Q (1.1) 
ev - (A - zI) H u = (1.2) 

and 

u H v = Q. (1.3) 

Following Theorem 1 1 . 1 1 and Lippert and Edelman 18 , Sections 4 and 5] we make the 
following assumption. 

Assumption 1.2. Assume A — zl satisfies the conditions of Theorem \1.1\ and 
that B = A — euv 11 is well- conditioned. That is, with z = a + f3i, the 2x2 matrix 

£aa Sap ^ s we n_ cona \itioned, where £ aa denotes the second partial derivative of 
£a/3 £pp J 

e with respect to a, etc. (see Theorem 5.1 and Corollary 5.2]). 

The paper is organised as follows. Section [2] contains some background theory 
and the derivation of the implicit determinant method for our problem. Section [3] 
describes the Newton method applied to this problem and in Section 13.11 we give 
numerical examples that illustrate the power of our approach. 
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2. The implicit determinant method to find a nearby defective ma- 
trix. In this section we describe some background theory and present our numerical 
approach to finding a nearby defective matrix, which is formulated as solving a 3- 
dimensional real nonlinear system. We emphasise that, since our numerical method 
uses standard Newton's method to solve the nonlinear system, we cannot guarantee to 
find the nearest defective matrix. However, a more sophisticated nonlinear solver may 
be used if greater reliability were sought. We do not do this here because in all our 
examples the nearest defective matrix was found using standard Newton's method. 

First, we formulate the problem following Alam and Bora [TJ Section 4] . Equations 
(|1.1[) - (|1.2[1 can be written as 



-el A-zI ' 




u 




' " 


(A - zI) H -el 




V 








(2.1) 



Set z 



i/3, x 



u 

v 



and introduce the Hermitian matrix 



K(a,j3,e) 



-el 

{A-{a + ip)I) H 



A - (a + 0)1 
-el 



(2.2) 



Then dimkerK(a,/3,e) 



Clearly, x is both a right and left null vector of K(a,/3,e). The following Lemma 
follows immediately from Assumption 11.21 

Lemma 2.1. Let e > satisfy the conditions in Theorem Furthermore, let 

u 

z = a + if3 be such that K(a, /3, e)x — 0, where x — 
1. 

We now introduce an algorithm to find the critical values of a, /? and e such that 
the Hermitian matrix K(a, /3, e) is singular and the constraint on x given by (|1.3|) is 
satisfied. We use the implicit determinant method, introduced in [9] to find photonic 
band structure in periodic materials such as photonic crystals. In [5] the implicit 
determinant method was used to find a 2-dimensional Jordan block in a Hamiltonian 
matrix in order to calculate the distance to instability. Here, in contrast to [S], we 
have a three-parameter problem with a constraint to satisfy, and apply the method 
to a classic problem in numerical linear algebra. 

First we introduce a bordered matrix M . The next theorem gives conditions to 
ensure that this matrix is nonsingular. 

Theorem 2.2. Let (a* , fi* ,e* ,x*) solve 

K(a t p,s)x = 0, x^O, 

so that dimkerK(a*,P*,e*) = 1 and x* 6 ker(K(a*, /3*, £*)) \ {0}. For some c £ C 2n 
assume 

c H x* j4 0. 



Then the Hermitian matrix 



M(a,P,e) 



K{a,0,e) 



(2.3) 



is nonsingular at a = a* , /? = /?*, e — e* . 

Proof. This result follows from [7J Lemma 2.8]. □ 
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As M(a*, (3* , £*) is nonsingular we have that M(a, f3, e) is nonsingular for a, (3 and e 
in the vicinity of a*, f3* and s*. Now consider the following linear system 



(2.4) 



where K(a,f3,e) is given by (|2.2I) . Clearly, Theorem 12.21 implies that both x and / 
are smooth functions of a, (3 and e near (a* , /?*,£*). 
Applying Cramer's rule to (|2.4[) we obtain 



K(a,j3,E) c 




x(a, /3, e) 




' " 


c H 




. /(a,/3,e) _ 




1 



f(a,(3,e) 



det K(a, (3, e) 
det Af (a, /3, e) 



and as M(a,/3,e) is nonsingular in the neighbourhood of («*,/?*,£*) by Theorem 
2.21 there is an equivalence between the zero eigenvalues of K(a,/3,e) (which we are 
looking for) and the zeros of f(a,f3,s). Hence, to find the values of a, (3 and e such 
that det K(a, f3, e) = we seek the solutions of 



f(a,/3,e) =0. 
If /(a*, /?*,£*) = 0, the first row of system (|2.4[) gives 

K(a*,/3*,e*)a;(a*,/3*,e*) = 0, 



(2.5) 
(2.6) 



that is, x(a* , (3* , e*) = x* is an eigenvector of if (a* ,/?*,£*) belonging to the eigenvalue 
zero. For the following derivation we use the notation 



x(a,/3,e) 



u(a,P,e) 
v(a,f3,e) 



(2.7) 



Note also that since K(a,f3,e) and M(a,(3,e) are Hermitian, f(a,(3,e) is real. Dif- 
ferentiating the linear system (12.41) with respect to a leads to 



K(a,/3,e) c 
c H 



x a (a,P,e) 
f a (a,f3,s) 



and the first row gives 



v(a,(3,e) 
u(a, /3, e) 




v(a,(3,s) 
u(a, (3, e) 



(2.8) 



(2.9) 



Multiplying this equation evaluated at (a* , (3* ,e*) from the left by the eigenvector 
x* H of K (a*, /3*, s*) gives 



f a (a*,p*,e*)= [ 



v 
it* 



= u u +w u = 2Ke(u w ), 



where we have used x* H c = 1 from (|2.4p . Similarly differentiating the linear system 
(|2.4|) with respect to (3 gives 



K(a, ,3, e) c 
c ff 



fp(a,P,e) 



v(a,(3,e) 
—u(a, (3, e) 




(2.10) 
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fp(a*,0*,e*)=i[ u* H v* H ] J u , =i(u* H v* 
Clearly 

f a (a*,/3*,s*) = and fp(a*, p*, e*) = 

Hence, we have reduced the problem of finding a solution to det K(a* , (3* , £*) = 
with u* H v* = 0, to that of solving g(a, /?, e) — 0, where 



V u ) 



-2lm(u* H v* 



u v — U. 



f(ot,/3,s) 
f a (a,f3,e) 
fp(a,/3,e) 



(2.11) 



which is three real nonlinear equations in three real unknowns. In the next section 
we describe the solution procedure using Newton's method. 

3. Newton's method applied to g(a,f3,e) — 0. In this section we describe 
how to implement Newton's method for the nonlinear system g(a,/3,e) = 0. We also 
obtain a nondegeneracy condition that ensures nonsingularity of the Jacobian matrix 
of g at the root, and hence confirms that Newton's method converges quadratically for 
a close enough starting guess. The nondegeneracy condition is shown to be equivalent 
to one introduced by Lippert and Edelman [5] for the conditioning of the 2-dimensional 
Jordan block of B = A — euv H . 

Newton's method applied to g(a,f3,s) is given by 



G(a w ,/3 W ,£ W ) 



7((a«j8« e«), 



(3.1) 



where o/ i+1 > = + Aa«, /3( 4+1 ) = /3« + A/3« and £( l+1 ) = £« + AeW, for i = 
0, 1, 2 . . . until convergence, with a starting guess (a^-°\ (3^°\e^), where the Jacobian 
is 



(3.2) 



and all the matrix entries are evaluated at (a^,ft (l \e^). The values of /«, f [ a l) and 
/i are found using (|2.4j) . (|2.8p and (|2 . 10|) . For the remaining values we differentiate 
dl3I), dHU) and ([2~TU|) with respect to £, that is, 





r fW 


Ai) 

h 


Ai) - 


G(a w ,/? w ,£ w ) = 


Ai) 
J aa 


Ai) 
J aft 


Ai) 
J ae 




Ai) 
. J Pa 


Ai) 

J pp 


Ai) 
J Pe J 



K(a,P,e) c 




x £ (a, ft, e) 




x(a, /3, e) 


c H 












(3.3) 



K(a,f3,e) c 
c H 



f a£ (a,(3,e) 



v E (a,(3,e) +u a (a,(3,e) 
u e {a.,f},e) +v a (a,f3,e) 




(3.4) 
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and 



K(a,/3,e) c 
c H 



iv £ (a, /3, e) + up(a, (3, e) 
-iu e (a,{3,e) +vp(a,P,e 




in order to find fg , fas and fj^. Furthermore, we differentiate 
a and /3, that is 



K(a,(3,e) c 
c H 



%aa (^5 ^) 

f aa (a,p,e) 



v a (a,(3,e) 
u a (a,/3,e) 




(3.5) 
with respect to 

(3.6) 



and 



K{a,p,e) 



x a p(a,f3,e) 
f a p(a,/3,e) 



iv a (a,(3,e) +vp(a,P,e) 
-iu a (a, (3, e) + up(a, /3, e) 




(3.7) 



to compute faa and j a L = /V. Finally, differentiate (|2.10|) with respect to /3 to get 



K(a,f3,e) c 
c H 



f/3p(a,f3,e) 



2 1 



-u/3(a,(3,e) 




(3.8) 



Therefore, in order to evaluate the components of G(a^ , /?W ; gM) anc [ g(a;W , /?W , gW) 
we only need to solve the linear systems above, which, importantly, all have the 
same Hermitian system matrix M(ofi- l \ I3^\e^). Hence only one LU factorisation of 
per iteration in Newton's method is required. Note that Newton's 
method itself is only carried out in three dimensions. Next, we state the Newton 
method algorithm for this problem. 

Algorithm 3.f (Newton's method). Given (a (0) , /3 (0) , e (0) ) and c £ C" such 
that M (a(°\ (3(°\ e^) is nonsingular; set i = 0: 
(i) Solve \2.$ and \2. 8\) and V2.10\) in order to evaluate 



/ a (aW (/ 8W,eW) 
/„ (a W, /?«,£«) 



(m,) 5ofce Iff. 3)) , fg. glh j3.7\ ), i3.8\) . \3.4\j and m order to evaluate the Jacobian 

G(aW,j8W,eW) gwen 6y EOI) . 
fjM^ Newton update: Solve $3.1]) in order to get (a^ l+1 \/3^ l+1 \s^ l+1 ^) 
(iv) Repeat until convergence. 

Finally we show, that provided a certain nondegeneracy condition holds, the 
Jacobian G is nonsingular at the root. In the limit we have 



G(a*, /?*,£*) 



f* e 

J aa J a/3 •> ae 

/# j** j** 

/3a J/3/3 J/3e 



(3.9) 



since /* = ft 
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Multiplying the first row of (I3.3j) evaluated at (a*,j3*,e*) from the left by x* H 



gives 



f s (a*,/3*,e*) = x* H x* > 0, since x* ± 0, 
(recall x* H c — 1 from (|2.4p V Hence the Jacobian (I3.9P is nonsingular if and only if 

(3.10) 



r a/3 ■— JuaJpfi ~ Ja/3 T u - 

With similar calculations as before we obtain 



f aa (a*,/3*,£*)=2x' 



* 



and 



.H 






< 


+ 




) 













(3.11) 



(3.12) 



Lemma 3.2. Under Assumption\TM F*p = ftJh ~ ftp* ^ °- 

Proof. If e is a simple singular value of (A — (a + f3i)I), qJgR, so that 

(A - (a + ,3z)7> = eit, (A - (a + /3i)I) H u = ev, 

then (see Sun 10J) e, u and v are smooth functions of a and f3. Furthermore, Lippert 
and Edelman Theorem 3.1] show that if u* H v* = then e* := e Q (a*,/3*) = 0, 
e*p := £,3(a*,/3*) = and B — A — eu*v* H has a 2-dimensional Jordan block. In 
addition, the ill-conditioning of the matrix B is determined by the ill-conditioning of 

see [51 Corollary 5.2]. Under Assumption 1 1 . 2 1 we have det(-B) ^ 0. 



E 



F F 

'-act °a/3 

_» 

e a/3 e pp 



Recall (|2.1[) and (|2.2I) where e = e(a,/3), u = u(a,/?), it = u(a,/3) and x 



u 
v 



Taking the second derivatives with respect to a and /3 and evaluating them at the 
root so that e*(a,/3) = e*p(a,(3) = we obtain 



K(a*,/3*,e*)x* 



u* 



K(a* , (3* ,e*)x%a + 2i 



EppX , 



and 



K( a *,/3*,e*)x* af) 



Multiplying those three equations by the eigenvector x* H of K(a* , (3* ,e*) from the 
left we obtain that 

faa = -i x * Hx *) £ *aa> ffip = ~{x* B X*)e*^ and f* p = -(x* H X*)e* a/3 , 

where we have used (|3~TTj) and ([3~T2]) . Hence F* = f* a f^-f*,3 2 ^ since det(E) ^ 
0. □ 
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In summary, Lemma l3. 21 shows that when the defective matrix B = A — euv H is well- 
conditioned Algorithm 13.11 should exhibit quadratic convergence for a close enough 
starting guess. 

Remark 3.3. We note that z = a + f3i is a saddle point of f (a, ft) and hence 

= faufpp ~ fap 2 < 0- This property can in fact be checked and is observed in all 
the computational examples in Section \3.1\ (see last column of Tables l^TiS. 5)) . 

We would like to note that our algorithm depends on the starting guess and hence 
does not guarantee convergence to the nearest defective matrix but only to a nearby 
one. However, all the algorithms currently known in the literature only find nearby 
defective matrices (see, in particular the methods suggested in [5]). 

We would also like to point out some computational advantages of our method. 
Both the method in [2] and our method provide a Newton method for finding a saddle 
point of a (in [5] and [5]) or /. For our problem the derivatives of / are particularly 
easy and simple to calculate. For any derivative a system with the same bordered 
Hermitian matrix (|2.3[) has to be solved - and we can get any 1st, 2nd or higher 
order derivatives by solving with the same matrix. Hence one matrix factorisation 
with costs of usually |(2n + l) 3 ~ -^n 3 or less for sparse systems is sufficient. Other 
explicit methods for calculating first and second derivatives have been derived (see [5] 
and [5]), which usually require a full SVD to be carried out, costing 21rt 3 operations 
(see [5]). Hence for large problems the implicit determinant method is more efficient. 
We show an example in the next section. 

3.1. Numerical examples. We now illustrate the numerical performance of 
our method with two examples which are taken from pQ. As has been mentioned 
earlier, since our method is based on Newton's method it finds a nearby defective 
matrix. We cannot guarantee it finds the closest defective matrix (this will depend on 
the starting guesses we use). However, in all cases considered here our method found 
the nearest defective matrix according to Alam and Bora pQ. 

Example 3.4. Let A e C" xn be the Kahan matrix fTTf , which is given by 

1 — c — c — c —c 
s —sc —sc —sc 

9 9 9 

J{ = S z —s z c —s^c 

s n-l 

where s Tl_1 = 0.1 and s 2 + c 2 — 1. We consider this matrix for n = 6, 15,20. As 
initial guesses we choose f)^ = and = for n = 6, — 0.12 for n — 15 and 
— 0.115 for n = 20. Further — cr mill7 u^ ' = u m i„ and = u m in> where cr min 
is the minimum singular value of A with corresponding left and right singular vectors 
it m in andv m i n . is determined from \2. 7\ ) and c = x^°\ We stop the iteration once 

\\g{a (l \(i {i \e {i) )\\<T, where r = 10~ 14 . 

Table [5~T1 shows the results for n = 6. In this case the eigenvalues 1.5849 x 10 _1 
and 10 _1 coalesce at 1.2763 x 10 _1 for a value of e = 4.7049 x 10~ 4 . The last column 
of Table O shows the value of = f&ff p - ( given by (|3.10[) ) and we see 
that the final value F*p ^ at the root. The quadratic convergence rate is clearly 
observed. 



DISTANCE TO A NEARBY DEFECTIVE MATRIX 



9 



Table 3.1 
Results for Example \3.4\ n = 6. 



i 


ad) 




e W 


||g(a«, || 













9.9694e-03 






1 


1.36436-01 





1.2145e-02 


8.1049e-02 


3.9318e-01 


2 


1.3319e-01 





7.1339e-04 


3.9165e-02 


-1.0032e+00 


3 


1.2767c-01 





4.9351e-04 


4.3976e-03 


-4.5529e-01 


4 


1.2763e-01 





4.7049e-04 


8.2870e-05 


-4.3191e-01 


5 


1.2763e-01 





4.7049e-04 


4.7344e-08 


-4.3136e-01 


6 


1.2763c-01 





4.7049e-04 


5.3655e-15 


-4.3136e-01 



Table 3.2 
Results for Example\3.4\ n = 15. 



i 


a« 




e (0 


||fl(aW,i9W,sW)|| 







1.2000c-01 





4.7454e-04 






1 


1.2042e-01 





2.1767e-06 


3.9203e-03 


-6.1848e-03 


2 


1.3116e-01 





1.0065e-06 


5.6943e-05 


5.6071e-06 


3 


1.2833e-01 





4.9786e-07 


2.8915e-05 


-6.7015e-05 


4 


1.2865e-01 





4.4839e-07 


1.6066e-06 


-5.9016e-05 


5 


1.2865e-01 





4.4850e-07 


1.7737e-08 


-6.1975e-05 


6 


1.2865e-01 





4.4850e-07 


1.9014e-12 


-6.1957e-05 


7 


1.2865e-01 





4.4850c-07 


3.5480e-18 


-6.1957e-05 



Table [321 shows the results for n = 15. In this case the eigenvalues 1.1788 x 10 
and 1.3895 x 1CT 1 coalesce at 1.2865 x 1CT 1 for a value of e = 4.4850e - 07 x 10~ 7 . 

Table 3.3 
Results for Example \3.4\ n = 20 . 



i 




0(0 


E« 


||g(c*W, || 


b aB 





1.1500c-01 





1.3141e-04 






1 


1.1507c-01 





1.1315e-07 


1.2702e-03 


-7.9071c-04 


2 


1.2010e-01 





3.2008e-08 


3.4299e-06 


-5.8539e-09 


3 


1.1997c-01 





1.8878e-08 


2.8840e-07 


-4.3105e-07 


4 


1.2000e-01 





1.9049e-08 


2.2944e-08 


-4.6343e-07 


5 


1.2000&-01 





1.9049e-08 


7.3704e-13 


-4.6360e-07 


6 


1.2000g-01 





1.9049e-08 


2.1281c-17 


-4.6360c-07 



Table |3~31 shows the results for n — 20. In this case the eigenvalues 1.1288 x 10 1 
and 1.2743 x lO" 1 coalesce at 1.2 x lO" 1 for a value of e = 1.9049 x 10~ 8 . 

From the last columns in Tables I3.HI3.3I we see that the value of F^p becomes 
smaller the larger the size of the Kahan matrix. This means the matrix B{e) — 
A — euv H becomes increasingly ill-conditioned as n increases. We also observe a 
corresponding deterioration in the rate of convergence of Newton's method as the 
value of F^i becomes smaller, which is consistent with the theory. 

Example 3.5. Let A e C nxn be the Grcar matrix taken from the Matlab gallery 
A = gallery ( 'grcar ' ,n), where n — 6,20. The eigenvalues of A appear in complex 
conjugate pairs and hence in this case two pairs of complex eigenvalues of A coalesce 
at two boundary points of the pseudo spectrum. 

As initial guess for n = 6 we take — 0, fi^ = —1, = 0, = u m i n 
and = fmin, where it m j n and v m { n are left and right singular vectors of A — fi^il, 
corresponding to the smallest singular value, x^ is determined from j2.7\ ). The 
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stopping condition is the same as in Example \S.4\ For n — 20 we take ft® = -2.5, 
the initial guesses for the remaining values are determined similarly. Furthermore 
c = x(°l 

Table 3.4 
Results for Examvle \3 '. 51 n = 6. 



i 


a« 






||ff(aW (( 8(0 )e (*))|| 










-1.0000e+00 









1 


1.2141c+00 


-2.3756e+00 


7.4297c-01 


5.0533c-01 


1.4186e-01 


2 


1.1159e+00 


-1.4291e+00 


9.5425e-02 


2.2193c+01 


-2.7279e+04 


3 


1.0512e+00 


-1.9848e+00 


4.3767c-01 


5.2914e-01 


-5.0768e+00 


4 


8.0543c-01 


-1.5940e+00 


1.4858e-01 


4.1255e-01 


-1.1717c+00 


5 


7.5742c-01 


-1.5944e+00 


2.1279e-01 


8.6847e-02 


-1.1323e+00 


6 


7.5335c-01 


-1.5912e+00 


2.1516e-01 


5.5621e-03 


-9.7810e-01 


7 


7.5332e-01 


-1.5912e+00 


2.1519e-01 


4.2790e-05 


-9.6333e-01 


8 


7.5332e-01 


-1.5912e+00 


2.1519e-01 


2.4851e-09 


-9.6323e-01 


9 


7.5332c-01 


-1.5912o+00 


2.1519c-01 


1.5798o-16 


-9.6323e-01 



Table [3T4l shows the results for n — 6. The eigenvalue pairs 1.1391 ± 1.2303z 
and 3.5849 x lO" 1 ± 1.9501* coalesce at 7.5332 x 10" 1 ± 1.5912i for a value of e = 
2.1519 x 1CT 1 . 

Table 3.5 
Results for Example \3. 51 n = 20 . 



i 


a« 


0{i) 




||fl(aW,/9W,eW)|| 










-2.5000c+00 











1 


9.5854e-02 


-2.3299e+00 


1.7989e-02 


1.3806e-01 


9.9103e-01 


2 


1.3904e-01 


-2.2465e+00 


1.3564e-03 


3.2308e-02 


-2.3623e-01 


3 


1.6141c-01 


-2.2042e+00 


7.2914e-04 


1.1930e-02 


-1.5963e-01 


4 


1.5554e-01 


-2.1818e+00 


4.5435e-04 


3.4851e-03 


-2.7982e-02 


5 


1.5338e-01 


-2.1815e+00 


4.9060e-04 


3.4265e-04 


-2.4693e-02 


6 


1.5331c-01 


-2.1817e+00 


4.9141e-04 


2.3240e-05 


-2.3956e-02 


7 


1.5331c-01 


-2.1817e+00 


4.9141e-04 


1.6942e-08 


-2.4012e-02 


8 


1.5331c-01 


-2.1817e+00 


4.9141e-04 


4.6672e-14 


-2.4012e-02 


9 


1.5331c-01 


-2.1817c+00 


4.9141c-04 


4.5263c-17 


-2.4012e-02 



Table l3~5l shows the results for n = 20. The eigenvalue pairs 1.0802xl0 -1 ±2.2253i 
and 2.1882 x lO" 1 ± 2.1132? coalesce at 1.5331 x 10" 1 ± 2.1817i for a value of e = 
4.9141 x 10~ 4 . 

The last columns in Tables 13.4113.51 show the values of F^l which converge to 
values away from zero. The latter iterates illustrate almost quadratic convergence. 
Note that in this example /3 ^ 0, so z = a + f3i is complex, though this makes no 
difference to the numerical method. 

We finally give a comparison of the method in [5] (see also [5]) to our method in 
terms of computational cost. Note that both the implicit determinant method - as all 
other methods known so far - only compute a nearby defective matrix. 

Example 3.6. Consider an n x n matrix with n — 1000, which is an identity 
matrix apart from the upper left 6x6 block which is the Kahan matrix from Example 
\3.4\ As initial guess we take the estimate which was used in J1J/. 

Table 13.61 shows the results for this comparison. We see that both the method 
from [2] and our new method exhibit very fast quadratic convergence to the desired 
nearby (in this case nearest) defective matrix (cf Table IBTTj) . However, the CPU times 
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Table 3.6 

Results for Examvle \3.6l n = 1000, Implicit determinant method (left) and method in fS? (right). 



i 


E W 


«W 


\\g(\W,sW)\\ 





4.6081e-04 


1.3175e-01 




1 


4.8049e-04 


1.2753e-01 


2.3311c-03 


2 


4.7050e-04 


1.2763c-01 


5.9278e-05 


3 


4.7049e-04 


1.2763&-01 


5.7187e-08 


4 


4.7049e-04 


1.2763e-01 


4.4987e-14 



i 


E W 


«W 


|| fl (AW,eW)|| 





4.6081e-04 


1.3175e-01 


4.6623e-03 


1 


4.7049e-04 


1.2753e-01 


1.1568e-04 


2 


4.7049e-04 


1.2763e-01 


5.6904e-08 


3 


4.7049e-04 


1.2763e-01 


1.3769e-14 



are very different. Whereas the method in [2] (right Table) requires a CPU time of 
24.3s the Implicit Determinant method only needs 5.4s. 

In summary, we note that both the method in [2] and the method described in 
this paper do not guarantee convergence to the nearest defective matrix. For large 
problems the Implicit Determinant method seems to be faster, as it is not necessary 
to compute the full SVD at each step. 

We note that we also compared the method described in [2] with the method 
described in this paper for Example 13.51 For this problem it is particularly hard to 
find good starting values in order for both methods to converge. If we generate the 
starting guesses as described in [2] we found that both methods stagnate or diverge 
- as for a small singular value e the Hessian becomes increasingly ill-conditioned. 
If we start with the starting guess described in Example 13.51 we found the Implicit 
Determinant method to converge (see Tables 13.41 and 13. 5p but the method described 
in [2] not to be defined as the second derivatives of the singular values are undefined 
for e^ ' = 0. However, to give a fair comparison, the method described in [5] describes 
a variant of Newton's method for the computation of a nearby defective matrix that 
is applicable to both generic and non-generic cases, whereas the method described in 
this paper only deals with the generic case, that is the computed singular value is 
assumed to be simple. 

4. Final remarks. We have developed a new algorithm for computing a nearby 
defective matrix. Numerical examples show that this new technique performs well 
and gives quadratic convergence in the generic cases. 

Also, since the method is only Newton's method on a real 3-dimensional nonlinear 
system (with only one LU factorisation required at each step) it is simple to apply 
and is significantly faster than the technique in pQ. 

However, as has already been mentioned, since it is based on Newton's method, 
convergence to the nearest defective matrix cannot be guaranteed, though in fact, in 
all the examples considered, convergence to the nearest defective matrix was achieved. 
Of course, a more sophisticated nonlinear solver, e.g. global Newton's method or a 
global minimiser, could be applied to (|2.11[) if required. 

Though our algorithm is designed to compute a nearby defective matrix in the 
generic case (that is, there is a well-conditioned 2-dimensional Jordan block) it has 
two features that enable it to recognise when the conditions of Assumption 1.2 fail. 
First, if there is another singular value near e then the condition number of M will 
be large. Second, if the condition number of M is small, but F a p is close to zero 
at the root, then this indicates the presence of a nearby matrix with a Jordan block 
of dimension greater than 2. As such the algorithm in this paper could be used to 
provide starting values for an alternative algorithm that could detect a higher order 
singularity. 
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